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ABSTRACT 

In this work, we investigate how and to which extent a quantum system can be driven along a prescribed path 
in space by a suitably tailored laser pulse. The laser field is calculated with the help of quantum optimal control 
theory employing a time-dependent formulation for the control target. Within a two-dimensional (2D) model 
system we have successfully optimized laser fields for two distinct charge transfer processes. The resulting laser 
fields can be understood as a complicated interplay of different excitation and de-excitation processes in the 
quantum system. 
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1. INTRODUCTION 

Since the first realization of the laser by T.H. Maiman in 1960 physicists and chemists dream about coherently 
controlling quantum systems using laser fields. For example, laser pulses may be applied to create and break 
a particular bond in a molecule, to control charge transfer within molecules, or to optimize high harmonic 
generation. The first approach to break a certain bond in a molecule using a laser tuned on resonance and 
initiating a resonance catastrophe failed. The molecule was distributing the energy internally too quickly, so 
that the specific bond did not break but instead the whole molecule was "heated". To overcome this so-called 
internal vibrational relaxation (IVR) a smarter excitation strategy and further technological improvements were 
necessary. 

With the advent of femtosecond laser pulses in the 1980's and the sophisticated technology^ for shaping 
the laser pulses the goal of controlling a complex chemical reaction with coherent light was finally achieved: 
For example, in 1998 Assion et al? showed that the product ratio CpFeCOCl+/FeCl+ of the organo- metallic 
compound (CpFe(C0)2Cl) can be either maximized or minimized by a specially tailored light pulse; or in 2001, 
Levis et al.^ demonstrated a rearrangement of molecular fragments. In both of these experiments adaptive laser 
pulse shaping techniques''^ have been applied, i.e., a computer analyzes the outcome of the experiment and 
modifies the laser pulse shape to optimize the yield of a predefined reaction product. This process is repeated 
until the optimal laser pulse is found. The number of experiments based on this so-called closed-learning-loop 
(CLL) is growing constantly, see for example Refs. 5-11. Recently, the pulse shaping techniques have been 
extended to allow also for polarization shaping,'^ i.e., experimentalists can indepentendly shape polarization, 
amplitude, and phase. 

Similar progress in quantum control can also be observed on the theoretical side starting in the 1980's with 
rigouros statements about controllability, '^^^ faster algorithms to calculate optimal laser pulses, ^'^'^ and a 
growing number of investigations for different types of systems and applications.'^' ^^^^ A large part of these 
studies employ a mathematical technique called optimal control theory (OCT)^^^^ well known from engineering. 
The employed numerical schemes have been designed to reach a predefined target at the end of a finite time- 
interval. Recently, these schemes have been extended to deal with time-dependent targets, meaning that we can 
also control the path the quantum system takes to the desired target. Little is known about the controllability^^^' 
of such objectives. 

In this article we employ the time-dependent target formalism to optimize laser fields that drive a charge along 
two distinct routes in 2D quadruple well. We allow the laser to have two independent components which would 
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require full polarization shaping in the experimental realization. The article is organized in the following way: 
In Sect. [2] we explain the basics of quantum optimal control theory for time-dependent targets and the algorithm 
which we use to optimize the laser field. The developed technique is then applied to a model system which is 
discussed in Sect. [3l The control objectives, i.e., the different pathways and the results of the optimizations are 
shown in Sect. [H We conclude the article with a summary and a short outlook in Sect. [5l 



2. THEORY AND ALGORITHM 

Let us consider an electron in an external potential V^(r) under the infiuence of a laser field. The laser pulse 
is assumed to propagate in z direction and therefore has two polarizations, x and so that a test charge is 
accelerated only in the plane perpendicular to the z axis. Given an initial state ^(r, 0) = 0(r) the time evolution 
of the electron is described by the time-dependent Schrodinger equation with the laser field modeled in dipole 
approximation (length gauge) 

i^vl/(r,i) = ^*(r,t), (1) 

H = H^-fie{t), (2) 
H, = f + V, (3) 

(atomic units [a.u.] are used throughout: h = m = e = 47reQ = 1). Here, fi = {jl^^fiy) is the dipole operator, 
and e{t) = (e^(t), ^y{t)) is the time-dependent electric field. For an electron we simply have fi = — r. The kinetic 

energy operator is T = — 

2.1. Derivation of control equations 

Our goal is to control the time evolution of the electron by the external field in a way that the time averaged 
expectation value of the target operator 0{t) is maximized. Mathematically, the optimization goal corresponds 
to the maximization of the functional 

1 '•^ 



JM = T^j^dt {^{t)\0{m{t)), (4) 

where 0{t) is assumed to be positive-semidefinite and defined by 

d{t) = d^{t) ^2TS{t-T)d^, (5) 

so we can also include targets in our formulation that only depend on the final time 7^.18,24,42 ^ examples 
will be discussed at the end of this section. 

The functional [^] will be maximized subject to a number of physical constraints. The idea is to cast also 
these constraints in a suitable functional form and then calculate the total variation. Subsequently, we set the 
total variation to zero and find a set of coupled partial differential equations. The solution of these equations 
will yield the desired laser field e{t). 

In more detail: Optimizing may possibly lead to fields with very high, or even infinite total intensity. 
In order to avoid these strong fields, we include an additional term in the functional which penalizes the total 
energy of the field. 



J^[e] = -a [ dt e^{t). 
Jo 



(6) 



The penalty factor a is a positive parameter used to weight this part of the functional against the other parts. The 
constraint that the electron's wave-function has to fulfill the time-dependent Schrodinger equation is expressed 

by 

J^[e,^,x] = -2Im dt (^xit)\{idt-H)\^it)). (7) 



with a Lagrange multiplier x(i*5^)- ^(1*5 ^) is the wave-function driven by the laser field e{t). The Lagrange 
functional has the form 



J[x,^,e]=J,m + J^[e] + J^[x,^,e]. (8) 
Setting the variations of the functional with respect to ^7 and e independently to zero yields 

ae^it) = -Im(x(t)|A,l*W), j = x,y,z (9) 

= {id,-H')^{r,t), (10) 

M/(r,0) = cPir), (11) 
{id,-H) x{r,t) + ^d,mir,t) = 

i(^x{r,t)-6^m{r,t))6{t-T). (12) 

Equation (|9|) determines the field from the wave-function ^(r, t) and the Lagrange multiplier x(r, t). While 
Equation (p!Q|) is a time-dependent Schrodinger equation for ^(r, t) starting from a given initial state 0(r) and 
driven by the field e{t). If we require the Lagrange multiplier x(r, t) to be continuous, we can solve the following 
two equations^^ instead of Equation (p!2|) : 

{id,-H)xir,t) = -^d,mir,t), (13) 
X(r,T) = d^^{r,T). (14) 

Hence, the Lagrange multiplier satisfies an inhomogeneous Schrodinger equation with an initial condition at 
t = T. Its solution can be formally written as 

X{r,t) = Ulx{r,to)-^ fdTUl(d,{T)^{r,T)), (15) 
where U^^ is the time-evolution operator defined as U^^ = Texp (^—i j^dt' H(t')^. 

The set of equations that we need to solve is now complete: (|9j), (p!Q|) . (pT]) . ([13]) and (p!4|) . To find an optimal 
field e{t) from these equations we use an iterative algorithm which is discussed in the next section. We conclude 
the derivation with a few examples for the target operator 0{t). 



Final-time control: Since the approach is a generalization of the traditional optimal control formulation given 
i]2i^'24,42 £^g|. observe that the latter is trivially recovered as a limiting case by setting 

6i(t)=o, 62 = -P = l*/)(*/l- (16) 

Here ) represents the target final state which the propagated wave-function ^(r, t) is supposed to reach at 
time T. In this case the target functional reduces to^^'^^ 

= (f (T)|P|vl/(T)) = |(vl>(r)|$^)|2. (17) 

The target operator may also be local in space. If we choose 0-y{t) = and O2 = S{r — Tq) (the density 
operator), we can maximize the probability density in Tq at t = T. 



J, = jdr {^{T)\d^\^{T))=n{r^,T). 



(18) 



Wave-function-follower: The most ambitious goal is to find the pulse that forces the system to follow a 
predefined wave- function <l>(r, t). If we choose 

d,{t) = \m){mi (19) 

62 = 0, (20) 

the maximization of the time-averaged expectation value of 0^{t) with respect to the field e{t) becomes almost 
equivalent to the inversion of the Schrodinger equation, i.e., for a given function <^>(r, t) we find the field e{t) 
so that the propagated wave- function ^(r, t) comes as close as possible to the target <^>(r,t) in the space of 
admissible control fields. 



Moving density: In this article we will focus on the following type of target. The operator used in Equa- 
tion (p!8|) can be generalized to 

d,{t) = S{v-r,{t)), (21) 



TJo 



= ^ I dt{^{mr-r,{tmit)) 



= ^ j^dtn{Y^{t),t). 

is maximal if the field is able to maximize the density along the given trajectory v^it). In Sect. [4] we will 
analyze the optimized fields calculated for two different trajectories of our 2D model system. In the numerical 
implementation the (5-function is approximated by a narrow Gaussian 

6(r) = 5{v-v^{t)) (22) 
1 r {v-voit)? ] 



■ exp 



2.2. Algorithm for time-dependent control targets 

Equipped with the control equations (|9]), (p!Q|) . (pT]) . (p!3|) . and (p!4|) . we have to find an algorithm to solve these 
equations for e{t). In the following we describe the scheme discussed in Refs. 37,40. A monotonic convergence 
in J can be proven if 77 G [0, 1] and ^ G [0, 2]?^ 

The algorithm starts with propagating ^^^^(0) = (j)^ forward in time with an initial guess for the laser field 

stepO: *(0)(0) ¥^){T). 

The backward propagation of x^^H^) started from x^^H^) — solving an inhomogeneous time-dependent 
Schrodinger equation which requires ^f^^\t) as input (with /c = 0), 

stepk: [<sm{T) (Iim(O)] . 

The brackets indicate that the storage of the wave function '^^^\t) can be avoided if we propagate it backwards 
in time as well using e^^\t). The backward propagation of x^^\t) requires the laser field determined by (/c = 0), 

^P{t) = (l-r,)ef (t)-^Im{x('=)(t)|A,|*('Hi)), J = ^.V- (25) 

The next step is to start a forward propag ation of ^(i)(0) = (j). [k = 0) 

[,(.)(0) 

[^(fe)(o) *('=)(r)] 

vEf(fc+i)(0) «'^^(*) *(fc+i)(T). 



and calculate the laser field 



ef+'\t) = {l-Oi'Ht)-l^lm{/'^^m^\¥'^+'Ht)), j = x,y. (26) 
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After the time evolution is complete we can close the loop and continue with Equation ([24|) . The whole scheme 
can be depicted symbolically by: 

stepO: ^(0)(0) '^'^ 
step k: 

(27) 

^(fe)(T)] 

The choice of and ^ completes the algorithm. ^ = 1 and r] = 1 correspond to the algorithm suggested in 
Ref. 42, while the choice ^ = 1 and 77 = is analogous to the method used in Ref. 18 with a direct feedback of 
^^^^(t). Further choices are discussed in Ref. 26. A more detailed discussion on the convergence of this algorithm 
with exclusively time-dependent target operators and a modified version of the target functional can be found in 
Refs. 41,43. An iteration scheme which also incorporates spectral restrictions^^ on the laser pulse is presented 
and analyzed in Ref. 41. 

Following the algorithm described above, one needs 5 propagations per iteration-step (avoiding the storage of 
the wave- function). Within the 2nd order split-operator scheme each time step requires 4 Fast Fourier Transforms 
(FFT),^^ because we have to know the wave function and the Lagrange multiplier in real space at every time-step 
to be able to evaluate the field from Eq. ([9]). This sums up to 2 * 10^ FFTs per 10^ time steps and iteration. 

3. MODEL AND COMPUTATIONAL DETAILS 

In this section we define our model system and analyze some of its properties. The potential surface is given by 

V{x, y) = a{x^ + y*) - b{x^ + y^) + c{x + y)\ (28) 

which describes a structure with a well in each corner of a square and a barrier in the center (see Fig. [2]). 
It is asymmetric with respect io x = —y^ i.e., the well at (—3,-3) is deeper than the one at (2.5,2.5), but 
symmetric for x = y. The parameters are chosen to be a = 1/64, 6 = 1/4 and c = 1/256 for which we show the 
six lowest eigenstates in Fig. [H In Table [1] we have calculated a few low-lying excitation energies necessary for 
the interpretation of the optimized pulse. 

4. RESULTS 

Finally, we employ the algorithm presented in Sect. [2] to calculate the laser field that guides the quantum system 
along the specified path. Both calculations start in the ground state |^(0) ) = |0q ) and employ the numerical 
parameters listed in Table O 

4.1. Transfer around the barrier 

The indirect route is described by 

i /2.5 - 5.5cos[7rV(2T)]\ 

i^oW-\^_3.0 + 5.5sin[7rt/(2T)]^ ' ^^^> 



The results of the optimization are presented in Fig. |31 Fig. |3(a)| shows that the x (in the upper panel) and 




Figure 1. (Color online). The six lowest eigenstates of the asymmetric quadruple well structure: (a) ground state, (b) 
1st, (c) 2nd, (d) 3rd, (e) 4th, (f) 5th excited state. The contour Unes indicate the potential surface. 



Table 1. Excitation energies in (a.u.) for the 2D asymmetric well, calculated by imaginary time propagation. 
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Table 2. Numerical parameters for the trajectory optimization of the two-dimensional quantum well structure. 
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Figure 2. The potential surface together with the "indirect" (a) and the "direct" (b) trajectory. The indirect route goes 
around the barrier while the direct route leads across the barrier. 



the y component of the field (in the lower panel) are quite different in their temporal behavior but similar 
in strength. The corresponding fluence [E^ = 1.0070) is rather high. The spectrum in Fig. |3(b)| contains 
significant contributions for Cc; G [0; 2] except for a narrow dip around = 0.5. This shows that a large number of 
intermediate states are occupied during the excitation process. Looking at the time evolution of the occupation 
numbers, we see that nearly all states up to the 10th excited state are significantly involved in the transition 
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Figure 3. (Color online). Results for the optimization of the indirect trajectory [see Fig. 2(a)| . (a) Optimized laser field, 
(b) Spectrum of optimized laser field. Lower panel: x component. Upper panel: y component, (c) Time evolution of the 

most important occupation numbers: | (n|^(t)) |^: Top panel: n = ( ), n = 4( ),n = 7 ( — • — ). Bottom panel: 

n = 5( ), n = 8 ( — • • — ), n = 9 (□). The ( ) line corresponds to the occupation for all states with n > 15. 

(d) Lower panel: Target trajectory [( ) Une], the expectation value of the position operator [(- - - -) line], and the 

position of the density maximum during the propagation (O ). Upper panel: Convergence of the algorithm where the 
(----) line corresponds to J, the ( ) line to Ji, and the laser fluence Eq is depicted by the ( ) Une. 



process. In Fig. 3(c)| we show the occupation numbers with the largest contributions, i.e., the ground state 



[( ) line], the 4th excited [(- - - -) line], the 5th excited [( ) line], the 7th excited [( — • — ) line], the 

8th excited [( — • • — ) line], and the 9th excited [(□) line]. We omit plotting the occupation number of the 3rd 
excited state since it is similar to the 4th. After t = 150 the pulse shifts nearly all occupation to the 5th exicted 
state which is localized at the target position in the right upper well. 

The convergence of the algorithm is shown in the upper panel of Fig. |3(d)[ We find a fast initial convergence 

(up to 20 iterations) and a slow improvement of the target functional [( ) line] for the remaining iterations. 

Although this scheme is in principle monotonically convergent it turns out to be very demanding to achieve an 
accurate enough solution of the control equations which would then guarantee the monotonicity. To assure the 
validity of the optimized laser field we have propagated it on a finer grid leading to the same results. 

The value of the functional after 200 iterations is = 0.0287 which appears to be small compared to the 
maximum possible value of 1. However, one has to be aware that is a measure for the shape of the density 
(which should be comparable to the target operator) and the position in time. So the control objective is quite 



demanding. To be able to better assess the quality of the pulse we compare the target trajectory with the 
expectation value of the position operator and the position of the density maximum, presented in the lower 
panel of Fig. |3(d)[ The figure shows that the expectation value meanders around the target curve rather closely. 
The position of the maximum is also close to the target trajectory but sometimes prefers to be in the minimum 
at (—2.7, 2.7) of the potential. Moreover, we analyze the time behavior of the density by plotting the two- 
dimensional density distribution for different points in time in Fig. [H Each figure contains a contour plot of the 
potential surface, the density distribution and a vertical line which marks the center of the target operator, i.e., 
the density is supposed to be maximized at this point. These figures illustrate that the laser guides the particle 
around the barrier in the desired way. 

4.2. Transfer over the barrier 

The trajectory leading directly across the barrier is given by 
with 0<t<T. 



The laser pulse shown in Fig. |5(a) is optimized to move the particle directly from the lower left well over the 
barrier to the upper right well. As one would expect, the x and y components of the laser pulse are identical 
in that case. Again the spectrum, shown in Fig. |5(b)[ is difficult to analyze. Compared to the spectrum for 



the indirect route [Fig. 5(a)| it is broader and therefore even more states are employed in the transition process. 



This can be seen in Fig. |5(c) , where we plot the occupation numbers for the most important states involved 



in the transition process. States up to the 20th excited state have a significant population. In contrast to the 
previous example the laser needs to excite many delocalized states to achieve a transition across the barrier. The 
convergence shown in the upper panel of Fig. |5(d)| is similar to the one found for the indirect route. After 200 
iterations we achieve a yield of = 0.0279 and a fluence of = 1.399 corresponding to a strong laser field. 
Again a deviation from the monotonic convergence is observed. 

In the lower panel of Fig. |5(d)] we plot the expectation value of the position operator [(- - - -) line] and the 
maximum of the density (O ) both indicate that the transfer of the particle occurs across the barrier. We can 
also see that a perfect localization cannot always be achieved by the laser, i.e., sometimes the position of the 
density maximum is found aside of the the target curve. 



5. SUMMARY 

We have demonstrated that it is possible to obtain laser pulses optimized to transfer a wavepacket along a 
predefined path. The laser fields have been calculated with the help of quantum optimal control theory using 
time-dependent targets. The pulses in our examples show that a complicated interplay of excitations and de- 
excitations is necessary to achieve the localization along the given trajectory. Both pulses involve a large number 
of eigenstates. Even for this rather simple system an optimization by hand using a combination of tt pulses^^ 
would be extremely cumbersome while the optimal control method requires only a naive guess for the initial 
laser field. 

In the example where the transfer occurs directly over the barrier the optimization converges to a linearly 
polarized pulse as one would expect from the symmetry of the problem. The shaping of linearly polarized 
laser fields is a rather established method. However, the optimized pulse for the transfer around the barrier 
requires a sophisticated time-dependent polarization. The experimental realization of polarization shaped pulses 
has been demonstrated recently. ^^"■'^^ Note that technological limitations can also be incorporated in the pulse 
optimization by combing the method presented above with additional restrictions on the optimal field, like for 
example spectral constraints.^^ 

Currently, we are working on the implementation of the optimal control algorithm into the freely available 
TDDFT-Solver package OCTOPUS^^'^'^ which in the future will allow us to investigate the controllability of 
multi-electron systems as well. 
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Figure 4. (Color online). Snapshots of the time evolution of the density in the potential well, (a) t = 0, (b) t = 40, (c) 
t = 80, (d) t = 120, (e) t = 160, (f) t = T = 200. The contour lines indicate the potential surface while the vertical line 
marks the center of the target operator [see Equation (|22)) ]. 




Figure 5. (Color online). Results for the optimization of the direct trajectory [see Fig. 2(b)| . (a) Optimized laser field. 

(b) Spectrum of optimized pulse. The x component [( ) line] and the y component [(- - - -) line] lie on top of each 

other, (c) Time evolution of the most important occupation numbers: | (n|^(t)) |^: top panel: n = ( ), n = 3 ( ), 

n = 9 (□), bottom panel: n = 5( ),n = 6 ( — • — ), n = 10 (O ). The ( ) line corresponds to the occupation 

for all states with n > 15. (d) Lower panel: (----) Expectation value of the position operator which lies on top of the 

target trajectory [( ) line], (O ) position of the density maximum during the propagation. Upper panel: Convergence 

of the algorithm. The (----) line corresponds to J, the ( ) line to Ji, and the ( ) Une to the laser fluence ^o- 
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